Causal mediation analysis is frequently used to assess potential causal mechanisms. More specifically, mediation analysis aims at disentangling the effects of a treatment on an outcome through alternative causal mechanisms and has become a popular practice in biomedical and social science applications. The aim of this project is to get familiar with mediation analysis and to understand its main sub-problems: indirect causal effect estimation and sensitivity analysis. You can either choose a more applied approach to apprehend this topic or choose the more theoretical paper on identification, inference and sensitivity.
Our work is based on the following articles:
LANGE T, VANSTEELANDT S, BEKAERT M. (2011). A Simple Unified Approach for Estimating Natural Direct and Indirect Effects. American Journal of Epidemiology. 176(3):190-195. https://academic.oup.com/aje/article/176/3/190/99496
IMAI K, KEELE L, YAMAMOTO T. (2010). Identification, Inference and Sensitivity Analysis for Causal Mediation Effects. Statistical Science. 25(1):51-71. https://www.jstor.org/stable/pdf/41058997.pdf?casa_token=ym7l4rAtIZkAAAAA:dDicuRrjqZW-dhh1c_aScoh8OMF7NtX5syaQXE82hYFImMQY_I3UB6DmuD8QpeBpaNRK5fYOxGtgjEad8CVBG4yfVyNzz8m87SYacIUgLOpSDa3HnJkNoQ
IMAI K, KEELE L, TINGLEY D, YAMAMOTO T. (2010). Causal Mediation Analysis Using R. Advances in Social Science Research Using R. VINOD H. D EDS. Vinod 129-154. https://cran.r-project.org/web/packages/mediation/vignettes/mediation-old.pdf
That Markown is complementray to the Latex report Causal Mediation Analysis and comes after it: we won’t repeat the bibliography for the sake of clarity, as well as the explanations of the methods alrdeady introduced in the report. The Markdown is divided into three parts:
library(mediation)
library(SparseM)
library(quantreg)
library(nlme)
library(mgcv)
library(foreign)
library(tidyverse)
library(dplyr)
library(rmarkdown)
library(missMDA)
library(stringr)
library(ggplot2)
library(gridExtra)
library(corrplot)
library(cdata)
library(shiny)
The dataset JOBS II comes from the Job Search Intervention Study (Vinokur and Schul, 1997). In the JOBS II field experiment, 1.801 unemployed workers received a pre-screening questionnaire and were then randomly assigned to treatment and control groups.
Those in the treatment group participated in job-skills workshops while those in the control condition received a booklet describing job-search tips:
Two key outcome variables were measured:
We will successively lead two analysis, (one for each outcome), one variable being continuous and the other discrete (binary).
In the JOBS II data, a key mediating variable is:
Finally, in addition to the outcome and mediators, the JOBS II data also include the following list of baseline covariates that were measured prior to the administration of the treatment:
data("jobs")
paged_table(jobs, options = list(rows.print = 10))
The study starts with depress2 for the outcome (continuous variable as well as the mediator job_seek).
We use the Barron-Kenny procedure to test mediational hypothesis in our causal model: * checking the correlation between the treatment and the outcome, * checking the correlation between the treatment and the mediator, (indeed, no mediation if no effect between the treatment and the mediator), * checking the impact of the correlation between the treatment and the mediator over the outcome.
As explained in our report, the first step is to estimate the regressions coefficients corresponding to the predictions of the mediator M jb_seek and of the outcome Y depress2, with sequential regressions (W referring to the treatment and X to the confounding factors).
\(Y=cW+\epsilon_1\)
\(M=aW+\epsilon_2\)
\(Y={c}'X+bW+\epsilon_3\).
Let’s estimate the coefficients.
model_1_y <- lm(depress2 ~ treat + depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income, data = jobs)
model_2_m <- lm(job_seek ~ treat + depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income, data = jobs)
model_3_y <- lm(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income, data = jobs)
Here is a summary of the last regression for the outcome.
summary(model_3_y)
##
## Call:
## lm(formula = depress2 ~ treat + job_seek + depress1 + econ_hard +
## sex + age + occp + marital + nonwhite + educ + income, data = jobs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.8247 -0.3834 -0.1090 0.2602 2.8287
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.4527281 0.1940639 7.486 1.74e-13 ***
## treat -0.0367886 0.0407940 -0.902 0.36740
## job_seek -0.1773802 0.0279535 -6.346 3.56e-10 ***
## depress1 0.4098612 0.0371003 11.047 < 2e-16 ***
## econ_hard 0.0679692 0.0221404 3.070 0.00221 **
## sex 0.0621566 0.0448206 1.387 0.16586
## age 0.0007858 0.0021942 0.358 0.72034
## occpmanegerial 0.0664879 0.0633665 1.049 0.29435
## occpclerical/kindred 0.0503458 0.0643692 0.782 0.43434
## occpsales workers -0.0348333 0.0836727 -0.416 0.67729
## occpcraftsmen/foremen/kindred -0.0290567 0.0800059 -0.363 0.71656
## occpoperatives/kindred wrks 0.1635053 0.0829922 1.970 0.04914 *
## occplaborers/service wrks -0.0215721 0.0848505 -0.254 0.79937
## maritalmarried -0.0072627 0.0551828 -0.132 0.89532
## maritalsepartd 0.2019970 0.1132861 1.783 0.07492 .
## maritaldivrcd -0.0453020 0.0639424 -0.708 0.47884
## maritalwidowed 0.0923133 0.1458613 0.633 0.52697
## nonwhitenon.white1 -0.1081444 0.0538549 -2.008 0.04494 *
## educhighsc -0.0023664 0.0901228 -0.026 0.97906
## educsomcol 0.0226457 0.0907839 0.249 0.80307
## educbach 0.0148269 0.1010784 0.147 0.88341
## educgradwk 0.1782504 0.1070053 1.666 0.09611 .
## income15t24k -0.0486597 0.0623912 -0.780 0.43565
## income25t39k -0.0208905 0.0641806 -0.325 0.74488
## income40t49k -0.0528838 0.0780279 -0.678 0.49811
## income50k+ -0.1179727 0.0735582 -1.604 0.10912
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.571 on 873 degrees of freedom
## Multiple R-squared: 0.2537, Adjusted R-squared: 0.2324
## F-statistic: 11.87 on 25 and 873 DF, p-value: < 2.2e-16
A quick glance at the p-values of the coefficients is enough to ensure that many null hypothesis of non-significance aren’t rejected at \(95\%\): age, marital status… It should also be noted that the the coefficient estimate of pre-level of depression is high (\(0.45\)).
However, we need to be careful in that analysis as the non-significance of a coefficient doesn’t mean that the variable has no impact on the outcome/mediator: that variable could be correlated to another variable concentrating the influence on the predicted variable. Furthermore, correlation does not imply causality AND no correlation does not disaprove causation (Bollen 1989), hence the causal mediation analysis (the causation can be in a mediation relation or hidden by confounding variables).
Here is a summary of the regression for the mediator.
summary(model_2_m)
##
## Call:
## lm(formula = job_seek ~ treat + depress1 + econ_hard + sex +
## age + occp + marital + nonwhite + educ + income, data = jobs)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.84605 -0.39919 0.05924 0.51764 1.46917
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.8806256 0.1947174 19.930 < 2e-16 ***
## treat 0.0774238 0.0492939 1.571 0.116624
## depress1 -0.2540256 0.0440638 -5.765 1.13e-08 ***
## econ_hard 0.1036040 0.0265612 3.901 0.000103 ***
## sex -0.0053180 0.0542355 -0.098 0.921913
## age 0.0005308 0.0026550 0.200 0.841575
## occpmanegerial 0.0056477 0.0766773 0.074 0.941302
## occpclerical/kindred -0.1132352 0.0777967 -1.456 0.145883
## occpsales workers -0.0137738 0.1012484 -0.136 0.891821
## occpcraftsmen/foremen/kindred -0.2015647 0.0965720 -2.087 0.037160 *
## occpoperatives/kindred wrks -0.2959024 0.0999259 -2.961 0.003147 **
## occplaborers/service wrks -0.3565544 0.1019639 -3.497 0.000494 ***
## maritalmarried 0.0381422 0.0667623 0.571 0.567934
## maritalsepartd 0.3641314 0.1365291 2.667 0.007793 **
## maritaldivrcd 0.2041101 0.0770659 2.649 0.008230 **
## maritalwidowed -0.3301324 0.1761481 -1.874 0.061240 .
## nonwhitenon.white1 0.0615794 0.0651346 0.945 0.344707
## educhighsc 0.1813264 0.1088818 1.665 0.096201 .
## educsomcol 0.1638371 0.1097146 1.493 0.135719
## educbach 0.2563072 0.1220038 2.101 0.035943 *
## educgradwk 0.2013935 0.1293041 1.558 0.119709
## income15t24k 0.1583888 0.0753071 2.103 0.035730 *
## income25t39k 0.0898314 0.0776033 1.158 0.247355
## income40t49k 0.1999402 0.0941763 2.123 0.034031 *
## income50k+ 0.1631108 0.0888391 1.836 0.066694 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.691 on 874 degrees of freedom
## Multiple R-squared: 0.1243, Adjusted R-squared: 0.1002
## F-statistic: 5.169 on 24 and 874 DF, p-value: 2.768e-14
It may seem interesting at first glance to look at the differences with the outcome: the coefficient for pre-level of depression (depress1) is significant and is negative (\(-0.25\)), instead of being positive when predicting the outcome (\(0.4\)), on the other side, variables previously declared as non-significant at \(95\%\) are now declared as significant (marital status). Nonetheless, contrary to basic correlation coefficients, regression coefficients can be affected by the coefficients of the other variables and we had better not rushing to draw conclusions.
Let’s to proceed to the mediation analysis instead. Actually, the package only requires in input the regression models estimated for the second and third equations (model_2_m and model_3_y).
out.1 <- mediate(model_2_m, model_3_y, sims = 1000, treat = "treat", mediator = "job_seek")
summary(out.1)
##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.0140 -0.0317 0.00 0.10
## ADE -0.0383 -0.1194 0.04 0.34
## Total Effect -0.0523 -0.1362 0.03 0.21
## Prop. Mediated 0.2210 -1.4844 1.82 0.27
##
## Sample Size Used: 899
##
##
## Simulations: 1000
The first column correspond to the point estimates and the next two columns to the estimated bounderies of confidence interval at \(95\%\). Here, the confidence intervals are estimated using quasi-Bayesian Monte Carlo approximations (\(1000\) draws). It makes sense as despite the fact that we provide fixed models to the function mediate, the estimation of the effects involve randomness (see equations (11)-(16) in the report).
If we look at the first three columns, it appears that the effects have a very low strength. One can easily confirm that with the the p-values, as they tend to indicate that none of the effects is significant.
Those results are reprensented in the graph below.
plot(out.1)
It appears that those effects (despite being non-significant) are more likely negative; if the job search mediates post-treatment depression, it’s negatively.
Eventually, one can check that - as expected - the Total Effect is the sum of the Average Diret Effect and of the Average Causal Mediation Effect (indirect effect), and that the proportion of the total effect due to mediation corresponds.
Even if that situation hasn’t been mentioned in the report, it could happen that the treatment moderates the causal mediation effect, i.e. the mediation effect could vary depending on the treatment status. The package mediation enables us to deal with that, using the Baron-Kenny procedure with Interaction term.
Actually, regarding the model it just comes down to consider a kind of mixed effect model, where the coefficients for the prediction of the outcome Y based on the mediator M change with the treatment W. (We can easily do that in R* as we are dealing here with linear models and a binary treatment)*. The model can be written as follows:
model_moder_y <- lm(depress2 ~ treat + job_seek + treat:job_seek +
depress1 + econ_hard + sex + age + occp +
marital + nonwhite + educ + income, data = jobs)
We display the results below. The intermediate effects have now two potential values, depending on whether the treatment is \(0\) or \(1\). Each time, both values are displayed, in dotted lines for the control group and in solid line for the treatment group.
out.2 <- mediate(model_2_m, model_moder_y, sims = 1000, treat = "treat", mediator = "job_seek")
summary(out.2)
##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME (control) -0.0189 -0.0448 0.00 0.10
## ACME (treated) -0.0120 -0.0281 0.00 0.10
## ADE (control) -0.0406 -0.1222 0.04 0.32
## ADE (treated) -0.0337 -0.1130 0.04 0.41
## Total Effect -0.0526 -0.1359 0.03 0.19
## Prop. Mediated (control) 0.2856 -1.6948 3.29 0.24
## Prop. Mediated (treated) 0.1875 -1.1292 2.24 0.24
## ACME (average) -0.0154 -0.0351 0.00 0.10
## ADE (average) -0.0371 -0.1153 0.04 0.35
## Prop. Mediated (average) 0.2365 -1.4329 2.91 0.24
##
## Sample Size Used: 899
##
##
## Simulations: 1000
plot(out.2, treatment = "both")
The effects don’t seem to be more significant in that framework than in the previous one. It happens that modelling a possible moderation in our situation doesn’t change anything, and it may not be necessary to go into that much details.
Considering linear models can sometimes be too restrictive. In that section we investigate options to study mediation effects with non-linear models. For that we consider Generalized Additive Models (GAMs) using splines to model the effects (piecewise polynomial functions).
model_gam_y <- gam(depress2 ~ treat + s(job_seek, bs = "cr") + depress1 + econ_hard +
sex + age + occp + marital + nonwhite + educ + income, data = jobs)
Besides, the quasi-Bayesian approximation is suited to linear models for the estimation of the confidence intervals (and more genrally parametric models) , but when it comes to non-linear models bootstrap happens to be the only reliable option (despite being more expensive in terms of computation ressources). We thus use non-parametric bootstrap to estimate confidence intervals in that section, resampling the data with replacement.
out.3 <- mediate(model_2_m, model_gam_y, sims = 1000, boot = TRUE, treat = "treat", mediator = "job_seek")
summary(out.3)
##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.0173 -0.0457 0.00 0.10
## ADE -0.0264 -0.0995 0.05 0.53
## Total Effect -0.0437 -0.1212 0.04 0.30
## Prop. Mediated 0.3953 -2.9543 3.78 0.35
##
## Sample Size Used: 899
##
##
## Simulations: 1000
plot(out.3)
Here no real added-value.
The mediator can also be modelled with a GAM of course, and it is also possible to deal with interaction terms using GAM’s.
Another option that hasn’t been mentioned in the report is to study the distribution of the outcome Y instead of its point values. More specifically, we look at the quantiles of that distribution.
Mediation modeling based on the mean may not be sufficient enough actually to understand different mediating effects across the outcome distribution. For that we replace the previous models for the outcome Y by models of quantile regression, using the package quantreg.
One should notice that the quantile \(tau\) has to be specify in advance. Here we consider the following values for \(tau: 0.2, 0.4, 0.6, 0.8\).
model_quant_y_2 <- rq(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income, tau = 0.2, data = jobs)
model_quant_y_4 <- rq(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income, tau = 0.4, data = jobs)
model_quant_y_6 <- rq(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income, tau = 0.6, data = jobs)
model_quant_y_8 <- rq(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age + occp + marital + nonwhite + educ + income, tau = 0.8, data = jobs)
out.4.2 <- mediate(model_2_m, model_quant_y_2, sims = 1000, boot = TRUE, treat = "treat", mediator = "job_seek")
out.4.4 <- mediate(model_2_m, model_quant_y_4, sims = 1000, boot = TRUE, treat = "treat", mediator = "job_seek")
out.4.6 <- mediate(model_2_m, model_quant_y_6, sims = 1000, boot = TRUE, treat = "treat", mediator = "job_seek")
out.4.8 <- mediate(model_2_m, model_quant_y_8, sims = 1000, boot = TRUE, treat = "treat", mediator = "job_seek")
Here is the summary of the effects for \(tau=0.2\).
summary(out.4.2)
##
## Causal Mediation Analysis
##
## Nonparametric Bootstrap Confidence Intervals with the Percentile Method
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.0121 -0.0283 0.00 0.084 .
## ADE 0.0419 -0.0494 0.08 0.550
## Total Effect 0.0298 -0.0590 0.07 0.838
## Prop. Mediated -0.4079 -6.5514 6.66 0.890
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 899
##
##
## Simulations: 1000
The plots are grouped below.
par(mfrow=c(2,2))
plot(out.4.2, main="Tau = 0.2")
plot(out.4.4, main="Tau = 0.4")
plot(out.4.6, main="Tau = 0.6")
plot(out.4.8, main="Tau = 0.8")
Only the quantile \(tau = 0.2\) differs from the others, the other results are quite the same as for the estimation with the mean. More specifically, the ACME doesn’t change and the only difference is with the ADE. Even if the ADE effect don“t appear to be significant and the x-axis of the plot \(tau = 0.2\) is more tightened than for the other plots (be careful with that), that effect is more positive than negative with \(tau = 0.2\).
One should notice that this result could hardly be attributed to a possible high variance of the effects for the tails of the distribution, provided that:
the ACME isn’t affected at all, only the ADE
the distribution on the left is quite flat, and the quantile \(tau=0.2\) may not be that instable
the other tail (way thiner) doesn’t face that issue.
The figure below may help support those arguments.
ggplot(jobs) + aes(depress2) + geom_density(alpha=0.6, color="darkolivegreen1", fill="darkolivegreen1") +
ggtitle("Symptoms of depression, post-treatment", subtitle="Density distribution")
We can thus affirm that - even if the effect isn’t significant - the treatment is more likely to increase depressive symptoms for the people almost not depressed at all than for the others. Nonetheless, one must keep in mind that those effects aren’t significant and that the confidence intervals don’t enable us to go further.
Apart from that quantile, we can nonetheless say that the mediating effects don’t vary much across the distribution of the outcome.
NB: we could also look at the distribution of the difference depress2 - depress1, symptoms of depression pre and post treatment (see below), as depress1 is taken into account in the regression.
ggplot(jobs) + aes(x=jobs$depress2 - jobs$depress1) +
geom_density(alpha=0.6, color="darkolivegreen1", fill="darkolivegreen1") +
scale_x_discrete(name="depress2 - depress1", limits=c(0:6)-3) +
ggtitle("Difference of symptoms of depression, pre and post treatment", subtitle="Density distribution")
Discrete outcomes can’t be treated as continuous outcomes (for the regressions amongs other things). Nonetheless, the package mediation handles that issue. We illustrate that with the discrete outcome work1.
To find out the estimate we use a probit regression instead of a linear regression. We could choose a logit regression, but the sensitivity analysis of the package mediation requires a probit regression.
model_discrete_y <- glm(work1 ~ treat + job_seek + depress1 + econ_hard + sex + age +
occp + marital + nonwhite + educ + income,
family = binomial(link = "probit"), data = jobs)
out.5 <- mediate(model_2_m, model_discrete_y, sims = 1000, treat = "treat", mediator = "job_seek")
summary(out.5)
##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME (control) 0.003149 -0.000933 0.01 0.146
## ACME (treated) 0.003387 -0.000982 0.01 0.146
## ADE (control) 0.055884 -0.001460 0.11 0.064 .
## ADE (treated) 0.056123 -0.001468 0.11 0.064 .
## Total Effect 0.059272 0.000880 0.11 0.050 *
## Prop. Mediated (control) 0.045554 -0.052445 0.44 0.184
## Prop. Mediated (treated) 0.049405 -0.053708 0.45 0.184
## ACME (average) 0.003268 -0.000964 0.01 0.146
## ADE (average) 0.056004 -0.001464 0.11 0.064 .
## Prop. Mediated (average) 0.047480 -0.052494 0.44 0.184
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 899
##
##
## Simulations: 1000
plot(out.5)
Although we have no interaction term in our model, we still have to a difference between the treatment and control groups, the gap can be overlooked as simply due to the non-linearity of the model.
Here , the idea of the sensitivity analysis is to study the robustness of our results to the ignorability assumption (and more precisely to the second condition as explained in the report). Indeed, conclusions that are harder to alter by a plausible violation of an identification assumption add a higher scientific value to analytics.
We keep considering \(1000\) simulations per call to the function mediate - as before - due to the variance of our results that require a high number of simulations (we have tried samller values but that wasn’t convincing enough).
Once again, we use the Baron-Kenny Procedure.
model_2_m <- lm(job_seek ~ treat + depress1 + econ_hard + sex + age + occp +
marital + nonwhite + educ + income, data = jobs)
model_3_y <- lm(depress2 ~ treat + job_seek + depress1 + econ_hard + sex + age + occp +
marital + nonwhite + educ + income, data = jobs)
As in the report, we present two ways of leading the sensitivity analysis: through the study of the residuals correlation \(\rho\), or through the coefficients of determination \(R\). The mediation package can do both.
med_cont <- mediate(model_2_m, model_3_y, sims=1000, treat = "treat", mediator = "job_seek")
sens_cont <- medsens(med_cont, rho.by = 0.05)
The rho.by option neables us to choose the precision for the incrementation of the parameter rho in the sensitivity analysis. Having a rough incrementation goes faster but don’t allow us to estimate precisely the robustness of our results. We choose the same value of \(\rho\) as in the article.
The object sens_cont contains the estimated effects and their CI’s corresponding to each value of \(\rho\) and \(R\).
summary(sens_cont)
##
## Mediation Sensitivity Analysis for Average Causal Mediation Effect
##
## Sensitivity Region
##
## Rho ACME 95% CI Lower 95% CI Upper R^2_M*R^2_Y* R^2_M~R^2_Y~
## [1,] -0.95 0.1808 -0.0418 0.4034 0.9025 0.5898
## [2,] -0.90 0.1183 -0.0274 0.2640 0.8100 0.5293
## [3,] -0.85 0.0894 -0.0207 0.1996 0.7225 0.4722
## [4,] -0.80 0.0715 -0.0166 0.1597 0.6400 0.4182
## [5,] -0.75 0.0588 -0.0137 0.1312 0.5625 0.3676
## [6,] -0.70 0.0489 -0.0115 0.1093 0.4900 0.3202
## [7,] -0.65 0.0410 -0.0096 0.0916 0.4225 0.2761
## [8,] -0.60 0.0342 -0.0081 0.0766 0.3600 0.2353
## [9,] -0.55 0.0284 -0.0068 0.0636 0.3025 0.1977
## [10,] -0.50 0.0232 -0.0057 0.0520 0.2500 0.1634
## [11,] -0.45 0.0185 -0.0047 0.0416 0.2025 0.1323
## [12,] -0.40 0.0142 -0.0038 0.0321 0.1600 0.1046
## [13,] -0.35 0.0102 -0.0030 0.0233 0.1225 0.0801
## [14,] -0.30 0.0064 -0.0025 0.0153 0.0900 0.0588
## [15,] -0.25 0.0028 -0.0026 0.0082 0.0625 0.0408
## [16,] -0.20 -0.0007 -0.0049 0.0036 0.0400 0.0261
## [17,] -0.15 -0.0040 -0.0105 0.0025 0.0225 0.0147
## [18,] -0.10 -0.0073 -0.0172 0.0026 0.0100 0.0065
## [19,] -0.05 -0.0105 -0.0242 0.0031 0.0025 0.0016
## [20,] 0.00 -0.0137 -0.0312 0.0037 0.0000 0.0000
## [21,] 0.05 -0.0169 -0.0382 0.0043 0.0025 0.0016
## [22,] 0.10 -0.0202 -0.0453 0.0050 0.0100 0.0065
## [23,] 0.15 -0.0234 -0.0526 0.0057 0.0225 0.0147
## [24,] 0.20 -0.0268 -0.0600 0.0065 0.0400 0.0261
## [25,] 0.25 -0.0302 -0.0677 0.0072 0.0625 0.0408
## [26,] 0.30 -0.0338 -0.0757 0.0080 0.0900 0.0588
## [27,] 0.35 -0.0376 -0.0841 0.0089 0.1225 0.0801
## [28,] 0.40 -0.0416 -0.0931 0.0098 0.1600 0.1046
## [29,] 0.45 -0.0460 -0.1027 0.0108 0.2025 0.1323
## [30,] 0.50 -0.0507 -0.1132 0.0118 0.2500 0.1634
## [31,] 0.55 -0.0558 -0.1247 0.0130 0.3025 0.1977
## [32,] 0.60 -0.0617 -0.1378 0.0144 0.3600 0.2353
## [33,] 0.65 -0.0684 -0.1528 0.0159 0.4225 0.2761
## [34,] 0.70 -0.0764 -0.1706 0.0177 0.4900 0.3202
## [35,] 0.75 -0.0862 -0.1925 0.0200 0.5625 0.3676
## [36,] 0.80 -0.0990 -0.2209 0.0229 0.6400 0.4182
## [37,] 0.85 -0.1169 -0.2609 0.0271 0.7225 0.4722
## [38,] 0.90 -0.1458 -0.3252 0.0337 0.8100 0.5293
## [39,] 0.95 -0.2083 -0.4647 0.0482 0.9025 0.5898
##
## Rho at which ACME = 0: -0.2
## R^2_M*R^2_Y* at which ACME = 0: 0.04
## R^2_M~R^2_Y~ at which ACME = 0: 0.0261
plot(sens_cont, sens.par = "rho")
We recall that the model is identifiable supposing that we know the correlation \(\rho\) between \(\epsilon_2\) and \(\epsilon_3\) (see report and part 1.2.1 of the Markdown).
Let’s explain what is depicted on the graph.
The solid line curve represents the mediation effect ACME against different values of \(\rho\), with the \(95%\) confidence interval (based on the delta method).
The dashed horizontal line represents the estimated mediation effect under the sequential ignorability assumption (previous quantification of the effect). It thus correspond to the effect we get when considering \(\rho = 0\) (with our current model).
The two extremities of the graph correspond to \(\rho = +/- 1\), i.e. a total correlation of \(\epsilon_2\) and \(\epsilon_3\) and thus (see part 1.2.1 and report) a total correlation of \(Y_i(w', m)\) and \(M(w)\) conditional to \(W=w\) (for any \(w, w', m\)), hinting thus at a very high mediation effect.
Why an \(AMCE\) positive when \(\rho\) becomes low enough (closer to \(-1\)) ? Because it corresponds to an unobserved pre-treatment covariate(s) that makes \(Y\) decrease when it makes \(M\) increase. By adjusting with regards to that covariate we would get that \(Y\) increases more (or decreases less) when \(M\) increases (through the intermediate effect).
How to interpret it ?
The critical point is when the curve crosses the x-axis: the mediation effect nullifies. That point corresponds here to \(\rho = -0.20\) (see summary above). Moreover, we can see that the \(95\%\) confidence intervals always contain \(0\), whatever the value of \(\rho\), and that our conclusion of non-significance of the \(AMCE\) effect may hold with any value of \(\rho\).
If we had had a significant effect previously, we would have looked at the range of \(\rho\) on which the \(95\%\) confidence interval \(CI\) don’t contain \(0\) in order to determine whether the results are robust. If that condition of \(0 \notin CI\) would only be true for a small interval of \(\rho\) around \(0\), that would be disappointed, otherwise it would hint at the fact that our results resist to a violation of the second condition required for ignorability.
We look at the products of the coefficients \(R^{2*}_M R^{2*}_Y\), corresponding to the proportion of the variance of our outcome which is unexplained by our current predictors, and that would be explained by a non-identified pre-treatment covariate (see report).
By the way, we have \(\rho^2 = R_M^{2*}R_Y^{2*}\) (see report) and \(\rho = sign(\lambda_2 \lambda_3) R_M^* R_Y^*\) with \(sign(\lambda_2 \lambda_3)\) the sign of the product; positive if the correlation with the unobserved covariate is of the same sign for the outcome and the mediator (see report), negative otherwise.
As above, the idea is to look at the estimated value of the \(AMCE\) for different values of \(R^{2*}_M\) and \(R^{2*}_Y\) that we choose, “what if the unobserved covariate had that precise importance ?”. More specifically the plot represents level curves of the \(AMCE\), (when the product is set).
par(mfrow=c(1,2))
plot(sens_cont, sens.par = "R2", r.type = "residual", sign.prod = "negative")
plot(sens_cont, sens.par = "R2", r.type = "residual", sign.prod = "positive")
The level curve corresponding to \(ACME = 0\) is in bold (critical point in the plot with \(\rho\)).
It is important to understand that the shape of the curves will always be the same as they correspond to \(xy= c^{te}\), for different constant values. (There are two plots to study both signs, as the two plots may not give the same values of ACME).
The point is thus to look at the value of ACME on each level curve: does our result hold when the product increases ? i.e. when we consider a neglected pre-treatment covariate with more influence, do we still get that the effect is significant ? The signs for extrem values can be explained as with \(\rho\).
This is only an graphical visualization of the summary, helping us to visualize when the ACME wouldn’t be really significant anymore.
Here the effect was already non-significant, and even with different values of the product it seems to remain non-significant. Nonetheless, if we originally had a significant effect, let’s say negative, we would have looked for the value of the product where the effect isn’t strictly negative anymore. The summary above can help. That would have indicated us the proportion of variance explained by neglected confounders which is required to have a situation where the effect would actually be null. That would help us having more hindsight on the robustness of our results by determining whether they would resist a violation of the ignorability assumption.
The objective of that part is to conduct a mediation analysis on our new dataset. We will use the methods explain in part 1.
We will consider the treatment acetalv, but as binary treatment: the groups “Some minorities” and “Many minorities” will thus be grouped together.
dd <- df_imputed %>% mutate_all(as.numeric)
dd$acetalv[dd$acetalv==2] <- 1
The outcome we chose is gvrfgap, but we will test all the mediators: imbleco, imtcjob, imwbcrm and imueclt. We select some of the confounding factors introduced below.
We quantify the effect with linear models, as non-linear models don’t bring any change (not shown here but we’ve looked at it). We have also looked at interaction terms and it makes very little difference too.
# mediator imbleco
model.m.1 <- lm(imbleco ~ acetalv + gndr + yrbrn + uempla + stfeco + tvtot +
happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)
model.y.1 <- lm(gvrfgap ~ acetalv + imbleco + gndr + yrbrn + uempla + stfeco + tvtot +
happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)
out.1 <- mediate(model.m.1, model.y.1, sims = 1000, treat = "acetalv", mediator = "imbleco")
# mediator imtcjob
model.m.2 <- lm(imtcjob ~ acetalv + gndr + yrbrn + uempla + stfeco + tvtot +
happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)
model.y.2 <- lm(gvrfgap ~ acetalv + imtcjob + gndr + yrbrn + uempla + stfeco + tvtot +
happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)
out.2 <- mediate(model.m.2, model.y.2, sims = 1000, treat = "acetalv", mediator = "imtcjob")
# mediator imwbcrm
model.m.3 <- lm(imwbcrm ~ acetalv + gndr + yrbrn + uempla + stfeco + tvtot +
happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)
model.y.3 <- lm(gvrfgap ~ acetalv + imwbcrm + gndr + yrbrn + uempla + stfeco + tvtot +
happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)
out.3 <- mediate(model.m.3, model.y.3, sims = 1000, treat = "acetalv", mediator = "imwbcrm")
# mediator imueclt
model.m.4 <- lm(imueclt ~ acetalv + gndr + yrbrn + uempla + stfeco + tvtot +
happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)
model.y.4 <- lm(gvrfgap ~ acetalv + imueclt + gndr + yrbrn + uempla + stfeco + tvtot +
happy + rlgblg + blgetmg + ipeqopt + impsafe, data = dd)
out.4 <- mediate(model.m.4, model.y.4, sims = 1000, treat = "acetalv", mediator = "imueclt")
We plot the results of the four mediators together.
par(mfrow=c(2,2))
plot(out.1, main="imbleco")
plot(out.2, main="imtcjob")
plot(out.3, main="imwbcrm")
plot(out.4, main="imueclt")
When looking at the summaries below, we can see that the three effects are almost always significant at \(95\%\) (one can look at the p-values but also at the confidence intervals on the plots that don’t contain \(0\)).
The Total Effect is relatively stable and the proportion of mediated effect really depends on the mediator.
We have for:
The rest is explained by the Direct Effect.
Based on the proportion of mediated effect in the total effect, we can distinguish between two groups of mediators:
imbleco and imwbcrm on the left: low mediation
imtcjob and imueclt on the right: explain almost one third of the total effect.
Here is the summary for imbleco: Taxes and services: immigrants take out more than they put in or less.
summary(out.1)
##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.00935 -0.01496 0.00 <2e-16 ***
## ADE -0.08498 -0.10511 -0.06 <2e-16 ***
## Total Effect -0.09432 -0.11577 -0.07 <2e-16 ***
## Prop. Mediated 0.09827 0.04271 0.16 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 40185
##
##
## Simulations: 1000
Here is the summary for imtcjob: Immigrants take jobs away in country or create new jobs.
summary(out.2)
##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.0259 -0.0311 -0.02 <2e-16 ***
## ADE -0.0680 -0.0877 -0.05 <2e-16 ***
## Total Effect -0.0939 -0.1147 -0.07 <2e-16 ***
## Prop. Mediated 0.2770 0.2140 0.37 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 40185
##
##
## Simulations: 1000
Here is the summary for imwbcrm: Immigrants make country’s crime problems worse or better.
summary(out.3)
##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.00443 -0.00953 0.00 0.076 .
## ADE -0.09002 -0.11038 -0.07 <2e-16 ***
## Total Effect -0.09445 -0.11472 -0.07 <2e-16 ***
## Prop. Mediated 0.04610 -0.00477 0.10 0.076 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 40185
##
##
## Simulations: 1000
Here is the summary for imueclt: Country’s cultural life undermined or enriched by immigrants.
summary(out.4)
##
## Causal Mediation Analysis
##
## Quasi-Bayesian Confidence Intervals
##
## Estimate 95% CI Lower 95% CI Upper p-value
## ACME -0.0293 -0.0357 -0.02 <2e-16 ***
## ADE -0.0657 -0.0847 -0.05 <2e-16 ***
## Total Effect -0.0949 -0.1162 -0.07 <2e-16 ***
## Prop. Mediated 0.3094 0.2373 0.39 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Sample Size Used: 40185
##
##
## Simulations: 1000
We conduct the same kind of sensitivity analysis as with the dataset JOBS II in part 1, i.e looking at the residual coefficient \(\rho\) and at the coefficients of determination \(R^*\).
med.cont.1 <- mediate(model.m.1, model.y.1, sims=1000, treat = "acetalv", mediator = "imbleco")
med.cont.2 <- mediate(model.m.2, model.y.2, sims=1000, treat = "acetalv", mediator = "imtcjob")
med.cont.3 <- mediate(model.m.3, model.y.3, sims=1000, treat = "acetalv", mediator = "imwbcrm")
med.cont.4 <- mediate(model.m.4, model.y.4, sims=1000, treat = "acetalv", mediator = "imueclt")
sens.cont.1 <- medsens(med.cont.1, rho.by = 0.05)
sens.cont.2 <- medsens(med.cont.2, rho.by = 0.05)
sens.cont.3 <- medsens(med.cont.3, rho.by = 0.05)
sens.cont.4 <- medsens(med.cont.4, rho.by = 0.05)
The residuals correlation plots are grouped below:
par(mfrow=c(2,2))
plot(sens.cont.1, sens.par = "rho", sub="imbleco")
plot(sens.cont.2, sens.par = "rho", sub="imtcjob")
plot(sens.cont.3, sens.par = "rho", sub="imwbcrm")
plot(sens.cont.4, sens.par = "rho", sub="imueclt")
With no surprise we can see that the higher the proportion of mediated effect, the more robust the result: indeed, for the mediators on the right, the effect remains strictly negative at \(95\%\) for a pretty wide interval of values \(\rho\) (the confidence intervals are also very thin). The confidence intervals for the first mediator imbleco are less thin but the effect remains strictly negative at \(95\%\) on a wide interval. The real issue regarding robustness concerns the mediator imwbcrm: the confidence intervals are large and contain zero for a a wide range of \(\rho\) values.
Besides, we can see on the summaries (below) that the ACME nullifies at \(\rho \leq -0.25\) for all mediators.
Here is the summary for imbleco: Taxes and services: immigrants take out more than they put in or less.
summary(sens.cont.1)
##
## Mediation Sensitivity Analysis for Average Causal Mediation Effect
##
## Rho at which ACME = 0: -0.25
## R^2_M*R^2_Y* at which ACME = 0: 0.0625
## R^2_M~R^2_Y~ at which ACME = 0: 0.0494
Here is the summary for imtcjob: Immigrants take jobs away in country or create new jobs.
summary(sens.cont.2)
##
## Mediation Sensitivity Analysis for Average Causal Mediation Effect
##
## Sensitivity Region
##
## Rho ACME 95% CI Lower 95% CI Upper R^2_M*R^2_Y* R^2_M~R^2_Y~
## [1,] -0.25 9e-04 -1e-04 0.0019 0.0625 0.0479
##
## Rho at which ACME = 0: -0.25
## R^2_M*R^2_Y* at which ACME = 0: 0.0625
## R^2_M~R^2_Y~ at which ACME = 0: 0.0479
Here is the summary for imwbcrm: Immigrants make country’s crime problems worse or better.
summary(sens.cont.3)
##
## Mediation Sensitivity Analysis for Average Causal Mediation Effect
##
## Sensitivity Region
##
## Rho ACME 95% CI Lower 95% CI Upper R^2_M*R^2_Y* R^2_M~R^2_Y~
## [1,] -0.95 0.0519 -0.0053 0.1092 0.9025 0.7630
## [2,] -0.90 0.0338 -0.0034 0.0711 0.8100 0.6848
## [3,] -0.85 0.0255 -0.0026 0.0535 0.7225 0.6108
## [4,] -0.80 0.0203 -0.0021 0.0426 0.6400 0.5410
## [5,] -0.75 0.0166 -0.0017 0.0348 0.5625 0.4755
## [6,] -0.70 0.0137 -0.0014 0.0288 0.4900 0.4142
## [7,] -0.65 0.0114 -0.0012 0.0240 0.4225 0.3572
## [8,] -0.60 0.0095 -0.0010 0.0199 0.3600 0.3043
## [9,] -0.55 0.0078 -0.0008 0.0163 0.3025 0.2557
## [10,] -0.50 0.0063 -0.0006 0.0132 0.2500 0.2113
## [11,] -0.45 0.0049 -0.0005 0.0103 0.2025 0.1712
## [12,] -0.40 0.0036 -0.0004 0.0077 0.1600 0.1353
## [13,] -0.35 0.0025 -0.0003 0.0052 0.1225 0.1036
## [14,] -0.30 0.0014 -0.0002 0.0029 0.0900 0.0761
## [15,] -0.25 0.0003 -0.0001 0.0008 0.0625 0.0528
## [16,] -0.20 -0.0007 -0.0014 0.0001 0.0400 0.0338
## [17,] -0.15 -0.0016 -0.0034 0.0002 0.0225 0.0190
## [18,] -0.10 -0.0026 -0.0054 0.0003 0.0100 0.0085
## [19,] -0.05 -0.0035 -0.0074 0.0004 0.0025 0.0021
## [20,] 0.00 -0.0044 -0.0093 0.0005 0.0000 0.0000
## [21,] 0.05 -0.0054 -0.0113 0.0005 0.0025 0.0021
## [22,] 0.10 -0.0063 -0.0133 0.0006 0.0100 0.0085
## [23,] 0.15 -0.0073 -0.0153 0.0007 0.0225 0.0190
## [24,] 0.20 -0.0082 -0.0173 0.0008 0.0400 0.0338
## [25,] 0.25 -0.0092 -0.0194 0.0009 0.0625 0.0528
## [26,] 0.30 -0.0103 -0.0216 0.0010 0.0900 0.0761
## [27,] 0.35 -0.0114 -0.0239 0.0012 0.1225 0.1036
## [28,] 0.40 -0.0125 -0.0263 0.0013 0.1600 0.1353
## [29,] 0.45 -0.0138 -0.0290 0.0014 0.2025 0.1712
## [30,] 0.50 -0.0151 -0.0318 0.0015 0.2500 0.2113
## [31,] 0.55 -0.0167 -0.0350 0.0017 0.3025 0.2557
## [32,] 0.60 -0.0183 -0.0386 0.0019 0.3600 0.3043
## [33,] 0.65 -0.0203 -0.0427 0.0021 0.4225 0.3572
## [34,] 0.70 -0.0226 -0.0475 0.0023 0.4900 0.4142
## [35,] 0.75 -0.0255 -0.0535 0.0026 0.5625 0.4755
## [36,] 0.80 -0.0292 -0.0613 0.0030 0.6400 0.5410
## [37,] 0.85 -0.0344 -0.0722 0.0035 0.7225 0.6108
## [38,] 0.90 -0.0427 -0.0898 0.0043 0.8100 0.6848
## [39,] 0.95 -0.0608 -0.1278 0.0062 0.9025 0.7630
##
## Rho at which ACME = 0: -0.25
## R^2_M*R^2_Y* at which ACME = 0: 0.0625
## R^2_M~R^2_Y~ at which ACME = 0: 0.0528
Here is the summary for imueclt: Country’s cultural life undermined or enriched by immigrants.
summary(sens.cont.4)
##
## Mediation Sensitivity Analysis for Average Causal Mediation Effect
##
## Sensitivity Region
##
## Rho ACME 95% CI Lower 95% CI Upper R^2_M*R^2_Y* R^2_M~R^2_Y~
## [1,] -0.3 -3e-04 -0.0012 6e-04 0.09 0.0651
##
## Rho at which ACME = 0: -0.3
## R^2_M*R^2_Y* at which ACME = 0: 0.09
## R^2_M~R^2_Y~ at which ACME = 0: 0.0651
We display the usual graph with one mediator per column, and negative correlations at the top and positive ones at the bottom.
par(mfrow=c(2,4))
plot(sens.cont.1, sens.par = "R2", r.type = "residual", sign.prod = "negative")
plot(sens.cont.2, sens.par = "R2", r.type = "residual", sign.prod = "negative")
plot(sens.cont.3, sens.par = "R2", r.type = "residual", sign.prod = "negative")
plot(sens.cont.4, sens.par = "R2", r.type = "residual", sign.prod = "negative")
plot(sens.cont.1, sens.par = "R2", r.type = "residual", sign.prod = "positive", sub="imbleco")
plot(sens.cont.2, sens.par = "R2", r.type = "residual", sign.prod = "positive", sub="imtcjob")
plot(sens.cont.3, sens.par = "R2", r.type = "residual", sign.prod = "positive", sub="imwbcrm")
plot(sens.cont.4, sens.par = "R2", r.type = "residual", sign.prod = "positive", sub="imueclt")
As with the residual coefficient plots, we can see that the effect remains strictly negative for positive products, and that it requires a certain amount of unexplained variance (explained by an unobserved confounder) to get an effect null in reality.
Graphically, one could guess that the corresponding square product for a null effect is around \(R_Y^{2*} R_M^{2*} = 0.1 \times 0.6 = 0.25^2\).
For all mediators except imwbcrm, the ACME is then still guaranteed to be strictly negative as long as the proportion of unexplained variance which can be explained by an unobserved confounder is less than \(25\%\) (for a confounder whose correlations with the mediator and the outcomes are of opposite sign, otherwise it’s okay). We also use the \(\rho\) to state that.
Our conclusions eventually depend on the mediator:
imbleco: the proportion of mediated effect is around \(10\%\), the mediation is significant (all the effects are significant), and the result can be considered as robust to the violation of ignorability (see sensitivity analysis),
imtcjob: the proportion of mediated effect is around \(25-30\%\), the mediation is significant and the result is robust,
imwbcrm: the proportion of mediated effect is around \(5\%\), with a mediation effect significant at \(90\%\) but not at \(95\%\), and when looking at the sensitivity analysis we can state that the result isn’t robust,
imueclt: the proportion of mediated effect is around \(25-30\%\), the mediation is significant and the result is robust.
As we can see it by looking at the meaning of those mediations, the main point lies in the causal effect between the treatment and the mediator. One may argue that this relationship hasn’t been studied here, but that would be incorrect as the tests we consider take that relation into account (effect \(a\) with the notations of the report).
Finally, we haven’t included the unused mediators as confounding factors when we studied each mediator in particular (too correlated: that would make the analysis too much tricky), and the unexplained variance in each study could (might) actually be attributed to the unused mediators.
Eventually, even if the proportions of mediated effect can be viewed as low, part of those proportions could actually sum when considering the whole set of mediators.
The two previous points are linked to the issue of multiple mediators raised in the report.
Let’s have a look at the meanings of the mediations, just to explain the meaning of the results we discuss.
A total effect which is negative (as in our situation) suggests that having minorities living in same current area would make you more inclined to agree with the fact that the Government should be more generous judging applications for refugee status.
The mediation effects we get as results (if proven to be robust enough) suggest that having minorities living in the current area would make people more “immigration friendly” through the idea that:
immigrants would take less than they put in (imbleco),
immigrants would create new jobs more than they would take (imtcjob)
immigrants wouldn’t make country’s crime problems worse (imwbcrm)
country’s cultural life is enriched by immigrants (imueclt).
In Europe, being in contact with minorities due to a geographical proximity would favor certain positive ideas as listed above, having people more “immigration friendly” (for the results proved as robusts).